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ABSTRACT. It has been known for a long time that the equivariant 2+1 wave map into the 2-sphere blows up 
if the initial data are chosen appropriately. Here, we present numerical evidence for the stability of the blow-up 
phenomenon under explicit violations of equivariance. 



1. Introduction 

The work presented here, is a continuation of the investigation of the wave map system carried out by the 
same authors. In [5] we presented our results on the blow-up in the equivariant case which we obtained 
by evolution of a 2 + 1 code which did not enforce the equivariance. As we pointed out this is already an 
indication that the blow-up phenomenon is stable against perturbations of the size of the truncation error of 
the evolution algorithm used. In the present paper we will give numerical evidence that this is also true for 
situations where initial data are used for which equivariance is broken explicitly already on the level of the 
initial data. 

For this investigation we have used the same setup as in [5]. So we will be brief in the description and for 
the details refer the reader to that paper. 

1.1. The wave map system. We use an extrinsic formulation of the 2+1 wave map, so we study maps U 
from 2 + 1 -dimensional Minkowski space M 2+1 into the unit-sphere, embedded into the Euclidean space 

K 3 

U : M 2+1 — > S 2 -+ R 3 

The unit-sphere is described as usual as the zero-set of the polynomial 0(z) = (z 1 ) 2 + (z 2 ) 2 + (z 3 ) 2 — 1 = 
5abZ A z s — 1, where we denoted the Euclidean metric on R 3 by 8 A b- This implies the restriction 

(1) <$>(U) = U A U A -\=0. 
on the map U . 

The wave map equations are obtained from an action principle using the action 

(2) £/[U,dU]= [ (d a U A d a U A +X<p{U)) dtdxdy. 

Jm 2 +' 

Here, (t,x,y) are Cartesian coordinates on M 2+1 and A is a Lagrange multiplier used to implement the 
constraint (1). 

Extremising (2) with respect to U A and A leads to the Euler-Lagrange equations 

D g U A - 21 U A = 
U A U A -l=0, 
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where □ is the usual d'Alembert operator □ = d tt — d xx — d yy . It is possible to eliminate the Lagrange 
multiplier with the help of the constraint equation. However, we choose not do this because in our numerical 
algorithm we solve the constraint and determine the Lagrange multiplier at every time-step. 

For the sake of clarity, we relabel the component functions of U as follows: u := U 1 , v := U 2 and w := C7 3 . 
Then, we can write the wave map system in the form: 



(3) 
(4) 



ii — d xx u — dyyU — 2Xu = 
V — d xx v — dyyV — 2Xv = 
w — d xx w — dyy.w — 2X w = 
(m, v, w) = u 2 + v 2 + w 2 — 1 = . 



This system has two non-trivial static solutions Us 

2x 



(5) 



u s {x,y) 



v s (x,y) 



2y 



■, ws(x,y) = ± 



l-x 2 -y 2 



l+x z +y z l+x z +y l l+x z +y z 

They describe the inverse of the stereographic projection to the sphere from the north- resp. south pole. 



1.2. Blow-up dynamics. The key feature to investigate the blow-up of the 2+1 dimensional wave map 
system is the scaling invariance of the equations (holds in all dimensions) and the energy (only in 2 + 1 
dimensions). This means that the equations as well as the energy are invariant under the transformation 

{t,x,y) — > (st,sx,sy) with set. 

Due to the fact that the energy is also scaling invariant, the 2+1 dimensional case is called the energy 
critical case. 

The first numerical results on the blow-up of the equivariant system were obtained by Bizoh et. al. in [ ]. 
The following three observations were made: First, when the energy of the initial data is too large then 
a singularity will form. Later, Sterbenz and Tataru [ ] specified in more detail under what conditions the 
2 + 1 wave map with the 2-sphere as target, has non-singular solutions. 

Second, close to the blow-up it is possible to rescale the dynamic solution U of (3) so that it approximates 
(in an appropriate Sobolev space) the static solution (5): 

(6) ]imU(t,s(t)x,s(t)y) = U s (x,y) 

t / i 

where T is the blow-up time and s(t) is the so called scaling function. The blow-up respectively the 
singularity formation appears as a shrinking of the rescaled static solution (5). This result was proven by 
Struwe [8]. In the same article was also shown that the existence of a non-trivial static solution is necessary 
for singularity formation. 

The scaling function s(t ) can be used to detect how the singularity formation proceeds. This was used in 
| ] and | ] for the numerical investigation of the blow-up. Bizori et. al. stated two properties for the scaling 
function: s(t) > for t < T (there is no solution for t > T anymore) and s(t ) \ for t /T. Raphael and 
Rodnianski [6] as well as Ovchinnikov and Sigal [ ] presented a detailed work on the blow-up dynamics. 
In both articles, an analytic form for the scaling function s(t ) was obtained. Ovchinnikov and Sigal were 
able to reduce the number of free parameters and therefore could give a more precise description of the 
scaling function s(t ). 

The third observation was that towards the blow-up the local kinetic energy at the point of the singularity 
formation goes to zero and the local potential energy approaches the value An, which is the energy of the 
static solution (5). 
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Those results were confirmed by Isenberg and Liebling [3]. In our previous work [ ] we also observed the 
expected blow-up behaviour. In addition we showed that the blow-up is stable under perturbations with a 
magnitude of the truncation error of the numerical scheme. 

1.3. Numerical setup. The numerical method which will be used to solve the equations (3) and (4) is 
the same as given in [ ]. Therefore, we will only briefly outline the most important points here. We 
discretise the spatial derivatives in the action functional using fourth order centred finite differences. This 
yields a semi-discrete action for finitely many degrees of freedom. The Euler-Lagrange equations for this 
action gives Hamiltonian equations of motion which are symplectic by construction, i.e., they preserve 
the canonical symplectic form of classical mechanics. The constraint (4) results in as many holonomic 
constraints as there are points on the numerical grid. In this way we obtain a Hamiltonian system with 
holonomic constraints. Our time integration method takes advantage of these properties: we use the Rattle 
method [1], a symplectic integrator for Hamiltonian systems with holonomic constraints. 

As described in [ ] we use the unit square Q, = [0, 1] x [0, 1] as our domain of integration with homogeneous 
Neumann boundary conditions at the outer boundary, i.e., on the sides with x = 1 and y — I, respectively. 
The other sides we regard as lines of symmetry, thus effectively enlarging the domain of computation to 
the square [—1,1] X [—1,1]. We use equal resolutions in both directions, i.e., Ax = Ay = 1/(N — 1) on a 
grid with N grid points in each direction. 

1 .4. Initial data. To guarantee that the constraint equation (4) is satisfied, the initial data are chosen as 

u (Q,x,y) = sin($o(V, a)) cos((p (a)) 
v {0,x,y) = sin(i3o(r, a)) sin(<po(a)) 
wo{0 7 x.y) = cos (#o(>, a)), 

where r and a are polar coordinates in R 2 related to the Cartesian coordinates x and y in the usual way 

x = r cos(a), y — rsin(a). 

'Equivariance' means that rotations in the -plane in M 2+1 are mapped to rotations around the z-axis 
in R 3 and is reflected in this parametrisation of the initial data as the requirements that 

tfo(r,cr) = #o(r), (po(a)=kG with^eZ. 

In our choice of initial data we explicitly break the equivariance by making #o depend on a but we keep 
the reflection symmetry across the lines x = resp. y = discussed above. The function #o(r, d) is defined 
as follows 

Ag(r)h{a) for re[r u r 2 } 

otherwise 

4 jr-n) (r 2 -r) l 4 
(r 2 - r\ ) (r 2 - r\ ) _ 

ho(a) for a e [0,ob] 

1 for a G [ffo, f - Co] 
>(f -a) for ae[f-ob,f] 

where ho(o) is a monotonically increasing function with ho(0) = B and fto(ob) = 1. It is chosen in such a 
way that h is at least ^ 4 , see Fig. 1. The function h{o) describes the deviation from equivariance, which 
would correspond to h{a) = 1. The parameter B measures the strength of the deviation (B = 1 corresponds 
to the equivariant case). We choose (po((j) — a as in the equivariant case with homotopy index k= 1. 



#o(r,CT) = 
S(r) = 

hip) - 



4 



JORG FRAUENDIENER AND RALF PETER 




Figure 1 . The angular perturbation function h{o). 



The initial data for the velocities are chosen as 



u(0,x,y) — d r u — cos i9o(r,cr) #n(r, a) - 

r 

v(0,x,y) = d r v = cos i> (r, a) i%(r, a) - 

r 

w(0,x,y) = d r w = -sini9o(r,a)^o(r,a) 



where #g(r) denotes the derivative of $o(r) with respect to its argument. These initial data describe a ring- 
shaped bump in the xy-plane around the origin. The choice of the velocity initial data results in a shrinking 
of this ring towards the origin. After the function w(t,x,y) reached its minimum close to the origin, the 
wave packet expands again. 



2. Scaling function 



As described above the blow-up dynamics is captured by the scaling function s(t) between the solution and 
the approximated static solution. In [5] we described how we determine this function in the equivariant 
case: since w is an axisymmetric function w(t,r) in that case we determine its second derivative with 
respect to r at r — at every time t and find .(f) as the appropriate factor between this and the second radial 
derivative of the static solution. 

When equivariance is broken then w is no longer axisymmetric and we need to extract the scaling function 
from the full matrix of second derivatives, the Hessian, of w at the origin. The Hessian Hs(t,x,y) of the 
rescaled static solution 

1 _ (r(.x,y)/ s(t) ) 2 

w* s (t,x,y) := w s (^)/, W ) = 1 , //^ • 
at the origin is proportional to the identity matrix: 

fls(f,0,0)=— 4rl 2 . 

s z {t) 

If the blow-up dynamics in the non-equivariant case is similar to the equivariant case then close to blow-up 
the off-diagonal terms of the Hessian of the solution will be small compared to the diagonal terms and we 
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can extract the scaling function s(t) simply from the trace of the Hessian using 

trtf s (f,0,0) = — 

s (0 

Alternatively, we could find s(t) also by taking the determinant of the Hessian using 

detfls(f,0,0) = -^r. 

s (0 

Geometrically this just means that we take either the mean curvature or the Gauss curvature at the origin 
of the surface defined by the graph of w as the indicator for the scaling function. It turned out that there are 
no essential differences so we used the Gauss curvature throughout. 



3. Blow-up results 



In Subsection 1 .2 the dynamics of the wave packet was described. Now we want to give some more details 
about the behaviour of the solution for initial data with large energy. If the energy of the initial data is 
large enough one expects a singularity formation as it is presented in [ ] and [ ]. In our setup, we interpret 
a change in the behaviour of w(t ,x,y) at the origin as the appearance of the blow-up. We have specified 
the initial data so that there remains a residual symmetry, namely the reflection symmetry with respect to 
both coordinate axes. This symmetry has the consequence that the value of the function w should remain 
constant, i.e., with our initial data w(f,0,0) = 1 throughout the evolution. However, if the energy is too 
large the numerical solution suddenly changes w(f,0,0) = — 1. From a geometrical point of view the 
solution switches from the stereographic projection from the south pole to the projection from the north 
pole, i.e., it suddenly approximates the other static solution. The reason for this seems to be due to the 
numerical method we are using. The Rattle method for the time integration always forces the solution 
onto the constraint manifold. For large energies it becomes increasingly difficult for the iterative projection 
algorithm to find a solution. It seems it is somewhat easier for the system to flip to the other solution. 

Using the previously introduced methods, we are now able to analyse the blow-up dynamics and singularity 
formation in the non-equivariant case. This will be analogous to the equivariant procedure presented in [ ]. 
For these simulations the value of the deviation from the equivariance B is fixed and the amplitude A is 
increased towards the critical value. As the indicator for the presence of the blow-up singularity we take 
the above mentioned flip from one static solution to the other. 

The first step in the analysis of the blow-up is the determination of the critical amplitude A*. Figure 2 shows 
the scaling function for different values of the amplitude A, where the parameter B is fixed to B = 0.8 and 
the numerical resolution isN ~ 1281. The qualitative behaviour is the same as in the equivariant case. An 
increase of the amplitude towards the critical value A* w 0.87150780 leads to an increase in the time the 
system remains in the quasi-static state. The appearance of this quasi-static, hovering state is due to the 
limited spatial resolution. If the number of grid points for the simulations is increased, the critical behaviour 
moves to higher values of the amplitude and reduces the duration of the hovering state. Therefore, our 
calculations are limited by the spatial resolution. 

The next step is to determine the blow-up time T. This can be found by fitting the last sub critical scaling 
function to the analytic expression (see [ ]) 

(7) s(f ) = — (T - 1) exp (- V- ln(T -»)+*). 

This formula was derived for the equivariant case but if the blow-up dynamics is a stable phenomenon 
then close to blow-up this formula should apply to the non-equivariant case as well. We fit the curve for 
A = 0.87150779 to function (7). This results in T = 0.93485135 for the blow up time and b = -2. 1435346 
for the parameter which depends on the initial data. The residual error for this fit was 8.0327838 • 10~ 9 . 
Figure 3 shows the result of the fitting procedure. The fit interval t € [0.865,0.8816] was chosen as a 
compromise between being close enough to the blow-up time (lower time bound) and the time domain, 
where the scaling function numerically converges (upper time bound). 
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Figure 2. The scaling function s(t) for different parameter values A and B = 0.8. The 
value A* = 0.87150780 is the critical amplitude. Computed with N = 1281. 




Figure 3. Fit of the scaling function s(t) to the analytical expression, given in [ ]. The 
computed blow-up is T = 0.94094524. The value A = 0.87150779 is the last under- 
critical value shown, i.e., for which the solution does not change its behaviour at the 
origin. Computed with = 1281. 



We present now, how the rescaled dynamic solution approximates the equivariant harmonic map near the 
blow-up. Since equivariance is broken this process is not isotropic anymore but instead depends on the 
angle a. In figure 4 we show the graphs of the component w(t,x,y) taken along the x-axis and along the 
diagonal. As one can clearly see, as time progresses the profiles agree increasingly better. This being the 
case, we show in figure 5 the successive stages of the rescaled dynamical solution taken only along the 
x-axis. 

To measure the deviation from the equivariant case, the difference between the two time steps when each 
of the wave packets reach their minimum is used. Additionally, the difference in the minima itself is used. 
Table 1 shows the numerical results of the two rescaled wave packets. The time f m ; n is defined as the 
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time when the wave packet along the x-axis respectively the line y = x, reaches its minimum w m i n (f m i n ). 
Based on those results, the relative deviation between the respective values are computed. Additionally, the 
respective minima vv m ; n (0) at t = are shown. 

x-axis diagonal rel. deviation 

~t~ 0.9296875 0.929375 3.3624748- 10" 4 

Wmin(fmin) -0.94862286 -0.94867828 5.8418118 • 10 5 
Wmin(O) 0.76663899 0.64367501 0.19103426 
Table 1 . Comparison of the ingoing wave packet along the x-axis and the diagonal 



From the results in table 1 and the graphs we conclude that the difference between the profiles along the 
x-axis and y = x vanishes during the evolution. The deviation along the two lines becomes significantly 
smaller, ending in nearly rotational symmetry. 

4. Conclusions 

We presented here for the first time numerical indications for a blow-up and singularity formation in the 
non-equivariant 2+1 dimensional wave map system. The methods which were developed for the analysis 
of the equivariant case of this system were extended and generalised to the non-equivariant case. It was 
possible to show that using initial data close to equivariance can also to a blow-up behaviour lead in the 
non-equivariant case. 

However, our simulations also showed the need for higher numerical resolutions in the 2-dimensional 
non-equivariant case to get a deeper and more detailed view into the blow-up dynamics and singularity 
formation. This can be done with grid refinement techniques or by changing the spatial discretisation of 
the equations completely. The use of (pseudo) spectral methods could be an appropriate way. 
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FIGURE 4. Time evolution of the wave packet for the amplitude A = 0.87150779. The 
slice along the x-axis and along the diagonal y = x is shown. 
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FIGURE 5. Ingoing wave packet for A = 0.87150779 and B = 0.8 rescaled with the 
scaling function s(t) for various values of t along the x-axis. The rescaled functions 
w(t,s(t)r) approximate the static solution ws(r). The function w(t,r) reaches its mini- 
mum at t = 0.9296875. 



